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Abstract. Based on mesoscale lattice Boltzmann (LB) numerical simulations, we investigate the effects of viscoelas¬ 
ticity on the break-up of liquid threads in microfluidic cross-junctions, where droplets are formed by focusing a liquid 
thread of a dispersed (d) phase into another co-flowing continuous (c) immiscible phase. Working at small Capillary 
numbers, we investigate the effects of non-Newtonian phases in the transition from droplet formation at the cross¬ 
junction (DCJ) to droplet formation downstream of the cross-junction (DC) (Liu & Zhang, Phys. Fluids. 23, 082101 
(2011)). We will analyze cases with Droplet Viscoelasticity (DV), where viscoelastic properties are confined in the 
dispersed phase, as well as cases with Matrix Viscoelasticity (MV), where viscoelastic properties are confined in the 
continuous phase. Moderate fiow-rate ratios ^(1) of the two phases are considered in the present study. Over¬ 
all, we find that the effects are more pronounced with MV, where viscoelasticity is found to influence the break-up 
point of the threads, which moves closer to the cross-junction and stabilizes. This is attributed to an increase of the 
polymer feedback stress forming in the corner flows, where the side channels of the device meet the main channel. 
Quantitative predictions on the break-up point of the threads are provided as a function of the Deborah number, i.e., the 
dimensionless number measuring the importance of viscoelasticity with respect to Capillary forces. 

PACS. 47.50.Cd Non-Newtonian fiuid fiows Modeling - 47.11.St Multi-scale methods - 87.19.rh Fluid transport 
and rheology - 83.60.Rs Shear rate-dependent structure 

1 Introduction 

Microfluidic devices are important in the studies that require 
control over droplet size in small scale hydrodynamics [1-6]. 

Common droplet generators used are T-shaped geometries [7, 

8] and flow-focusing devices [9-11]. In T-shaped geometries, 
a dispersed (d) phase is injected perpendicularly into a con¬ 
tinuous (c) phase [7]. Forces are created by the cross-flowing 
continuous phase which periodically break off droplets; the 
flow-focusing geometry, instead, creates droplets by focusing 
a fluid thread into another fluid [9,10,12-14]. The operational 
regime of these devices is characterized by the Capillary num¬ 
ber Ca, which quantifles the importance of the viscous forces 
with respect to the surface tension forces, and the flow-rate 
ratio Q of the two fluids. Previous research has been mainly 
restricted to Newtonian fluids. However, there are many situ¬ 
ations of practical interest which result in considering a non- 
Newtonian behaviour [15-18]. The formation of viscoelastic 
droplets in Newtonian continuous phases was investigated in 
various flow-focusing geometries by Steinhaus et al. [16]. The 
effects of elasticity on filament thinning was studied by Ar- 
ratia et al [15,18] using dilute aqueous polymeric solutions 
with different molecular weights. In these studies it is observed 
that elasticity prolongs the processes of thinning of the fila¬ 
ment and significantly increases the interval required for break¬ 
up to complete. In a recent paper, Derzsi et al [11] presented 
an experimental study of the effects of elasticity in microflu¬ 


idic flow-focusing devices: the authors And that the elasticity of 
the focusing liquid facilitates the formation of smaller droplets 
and leads to transitions between various regimes at lower ra¬ 
tios of flow and at lower values of the Capillary numbers in 
comparison to the Newtonian case. Complementing these re¬ 
sults with systematic investigations by varying deformation 
rates and fluid constitutive parameters would be of great in¬ 
terest [7-10,19]. In a recent paper, Liu & Zhang [9,10] per¬ 
formed numerical simulations of three-dimensional microflu¬ 
idic cross-junctions at low Capillary numbers. A regime map 
was created to describe droplet formation at the cross-junction 
(DCJ), downstream of the cross-junction (DC), and stable par¬ 
allel flows (PF). The influence of flow-rate ratio. Capillary 
number, and channel geometry was then systematically stud¬ 
ied in the squeezing-dominated DCJ regime. 

In the present paper we study the impact of viscoelastic phases 
in the transition from DCJ to DC regimes. At fixed flow condi¬ 
tions, these effects will be shown to be particularly relevant for 
the case where viscoelastic properties are confined in the con¬ 
tinuous phase, due to the fact that the action of the co-flowing 
liquid is supplemented with the polymer feedback stresses that 
are well pronounced in correspondence of the comers where 
the side channels meet the main channel. Quantitative details 
on how the polymer feedback stresses change the transition 
from DCJ to DC regimes are provided. 

The paper is organized as follows: in Sec. 2 we will present 
the necessary mathematical background for the problem stud- 
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ied, showing the relevant equations that we integrate in both 
the continuous and dispersed phases. In Sec. 3 useful bench¬ 
marks for the shear and elongational rheology of the numerical 
model will be provided for the typical parameters used in our 
study. In Sec. 4 we will present the numerical results and char¬ 
acterize the effects of viscoelasticity on the various regimes. 
Conclusions will follow in Sec. 5. 


2 Macroscopic equations 

Numerical modeling of viscoelastic fluids often relies on cou¬ 
pling constitutive relations for the stress tensor, typically ob¬ 
tained via approximate representations of some underlying 
micro-mechanical model for the individual polymer molecules, 
with a Navier-Stokes (NS) description for the solvent. We will 
use a hybrid algorithm combining a multicomponent Lattice- 
Boltzmann (LB) model with Finite Differences schemes (FD): 
the solvent part of the model is obtained with LB models [20- 
23], which proved to be extremely valuable tools for the sim¬ 
ulation of droplet deformation problems [24-27], droplet dy¬ 
namics in open [28,29] and conflned [27,30,31] microfluidic 
geometries. LB is instrumental to solve the diffuse-interface 
hydrodynamic equations of a binary mixture of two compo¬ 
nents [32-37]: the resulting physical domain can be partitioned 
into different subdomains, each occupied by a “pure” fluid, 
with the interface between the two fluids described as a thin 
layer where the fluid properties change smoothly. The solute 
part of the model is based on a FENE-P constitutive model [38, 
39]. The EENE-P constitutive equations are solved with a ED 
scheme which is coupled with the solvent LB as described 
in [40,41]. Our numerical approach has been extensively vali¬ 
dated in our previous works [40,41], where we have provided 
evidence that the model is able to capture quantitatively rheo¬ 
logical properties of dilute suspensions as well as deformation 
and orientation of single droplets in conflned shear flows. The 
main essential features of the model are recalled in the Ap¬ 
pendix. 

Eor the geometry studied in this paper, we will provide quan¬ 
titative details on how the viscoelastic model parameters affect 
the break-up properties. We will analyze separately cases with 
Droplet Viscoelasticity (DV), where the viscoelastic properties 
are conflned in the dispersed (d) phase undergoing the break¬ 
up process, as well as cases with Matrix Viscoelasticity (MV), 
where the viscoelastic properties are conflned in the continuous 
(c) phase. In the MV case, the equations we solve in the con¬ 
tinuous phase are the Navier-Stokes (NS) equations coupled to 
the EENE-P constitutive equations 


V (r/,(Vuc + (Vucf)) + ^V • [f{rp)C]. 

d,C+ K • V)C = C ■ ( Vm,)+( • C 

_ ( 2 ) 


Here, Uc and Pc are the velocity and the dynamic viscosity 
of the continuous phase, respectively, pc is the solvent den¬ 
sity, Pc the solvent bulk pressure, and {'VucY the transpose 


of (Viic). As for the polymer details, r(p is the viscosity 
parameter for the EENE-P solute, Tp the polymer relaxation 
time, C the polymer-conformation tensor, I the identity tensor, 
f{rp) = {L^ — 3)/{L^ — rp) the EENE-P potential that ensures 
flnite extensibility, rp = ^Tr{C) and L is the maximum pos¬ 
sible extension of the polymers [42,43]. In the dispersed phase 
we just consider the NS equations 


Pd [dfUd + {ud • V)ud] = -VPd 

+v (T]rf(VMrf + (VMrf)^)) (3) 

where the different flelds have the same physical meaning but 
they refer to the dispersed phase. Immiscibility between the 
dispersed phase and the continuous phase is introduced using 
the so-called “Shan-Chen” model [40,44,45] which ensures 
phase separation with the formation of stable interfaces be¬ 
tween the two phases characterized by a positive surface ten¬ 
sion ( 7 . 

Eor the DV case, we consider the reversed case, where the 
EENE-P constitutive equations are integrated in the dispersed 
phase (i.e. (l)-(2) with c ^ d), while only the NS equations are 
considered in the continuous phase (i.e. (3) with d ^ c). 

As for the geometry used, we consider the simplest case where 
channels have a square cross-section H xH. The square cross- 
section is resolved with H x H = 50 x 50 grid points, while 
the main and side channels are resolved with = 1150 and 
Ly = 250 grid points, respectively. Besides the geometrical pa¬ 
rameters, the Newtonian problem is described by six parame¬ 
ters characterizing the flow and material properties of the flu¬ 
ids. These parameters are the mean speeds of the continuous 
and dispersed phases, Vc and vj, respectively; the viscosities 
of the two fluids Pc and of Eqs. (1) and (3), the interfacial 
tension CT, and the total density pc = ps = P (the same for the 
dispersed and continuous phases). We will assume perfect wet¬ 
ting for the continuous phase, while the dispersed fluid does not 
wet the walls. As usual with these kinds of systems [7,9,10], 
we choose the following groups: the Capillary number calcu¬ 
lated for the continuous phase. 


Ca = 



( 4 ) 


the Reynolds number Re = pvcH/(rijoj^c)^ the viscosity ratio 
A, and the flow-rate ratio 


Qd 
c ~ Q 


where and Qc = VcH^ are the flow-rates at the two 

inlets. Eor the flow regimes under consideration, the Reynolds 
number is small (Re ~ 0.01 — 0.1), and does not influence the 
droplet size, which leaves us with three governing parameters: 
Ca, A and Q. Notice that the total viscosity in the continuous 
phase T7 tot,c is either if]joT,c = flc + flp (for MV) or tJtqt^c = Pc 
(for DV). In the outlet, we impose pressure boundary condi¬ 
tions and use Neumann boundary conditions for the velocity 
held. A Dirichlet boundary condition is imposed at the inlets 
by specifying the pressure gradient that is compatible with the 
analytical solution of a Stokes flow in a square duct [25]. As for 
the polymer boundary conditions, we impose a Dirichlet type 
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boundary condition by linearly extrapolating the conformation 
tensor at the boundaries. 


3 Viscoelastic parameters and rheology of the 
numerical model 


The viscosity ratio of the two Newtonian phases / rjc [40,41] 

is tunable, a fact that allows to work with unitary viscosity ra¬ 
tio, defined in terms of the total (fiuid + polymer) shear viscos¬ 
ity ^ = rid/{ric ^rip) = 1 for MV and A = (r]j + rip)/ric = 
1 for DV. Consistently, we will compare the non-Newtonian 
simulations with the corresponding Newtonian case at A = 
/^c = 1 • The ratio between the polymer viscosity and the 
total viscosity is chosen as Tip/{ric4 + ^p) ~ 0.40. In order to 
report data with dimensionless numbers, we choose to quantify 
the degree of viscoelasticity with the Deborah number that we 

define as De = i^ 

mal stress difference which develops in the viscoelastic phase 
in presence of a homogeneous steady shear [39,42]. In the def¬ 
inition of the Deborah number, the viscosity is obviously in¬ 
dicated in the viscoelastic phase, either for MV or 

for DV. The shear rheology of the model can be quanti¬ 
tatively verified in the numerical simulations. There are indeed 
exact analytical results we can get by solving the constitutive 
equations for the hydrodynamical problem of steady shear fiow, 
Ux = fy, Uy = Uz = 0: both the polymer feedback stress and the 
first normal stress difference Ni for the FENE-P model [39,42] 
follow 


-/MC„ = 

sinh 1 larcsinh [ 2L 
I 3 I 4 



( 6 ) 


Ni 


^f(rp)(C^-Cyy) = 8^ 
Tp Tp 



X 


sinh^ 


-arcsinh 

3 




( 7 ) 


where A = yrp is the dimensionless shear. The validity of both 
Eqs. (6) and (7) is benchmarked in Eig. l(a)-(d): numerical sim¬ 
ulations have been carried out in three dimensional domains 
with // X // X // = 20 X 20 X 20 cells. Periodic boundary con¬ 
ditions are applied in the stream-flow (x) and in the transverse- 
fiow (z) directions while two walls are located at y = 0 and 
y = H. The linear shear fiow Ux = fy, Uy = Uz = 0 is imposed in 
the numerics by applying two opposite velocities in the stream- 
fiow direction (Ux{x,y = 0,z) = —Ux{x,y = H^z) = Uy^) at the 
upper (y = H) and lower wall (y = 0) with the bounce-back rule 
[46]. We next change the shear in the range 10“^ < lU^/H < 
10~^ Ibu (lattice Boltzmann units) and the polymer relaxation 
time in the range 10^ < Tp < 10^ Ibu at fixed rfp. The vari¬ 
ous quantities are made dimensionless with the viscosity rjp 
and the relaxation time Tp, and they are plotted as a function 


of the dimensionless shear A = Tpf. The values of the confor¬ 
mation tensor are taken when the simulation has reached the 
steady state. As we can see from the figures, all the numeri¬ 
cal results well agree with equation (6) and (7). In particular, 
both the shear stress and normal stress difference increase at 
large A to exhibit a sublinear behaviour (see figures 1(a) and 
1(c)) which directly relates to thinning effects in the dimen¬ 
sionless polymer shear viscosity, f{rp)Cxz/A, and first normal 
stress coefficient, Wi = f{rp){Cxx — Cyy)/A^ (see figures 1(b) 
and 1(d)). In the limit of Hookean dumbbells (Oldroyd-B limit, 
^ I) we can use the asymptotic expansion of the hyperbolic 
functions and we get N\ ^ 2Tpr{py^, so that 


De = 


^p Vp 

ric,d + riP ' 


( 8 ) 


Equation (8) shows that De is clearly dependent on the ratio 
between the polymer relaxation time Tp and the time Th defined 
as 


% 


H{ric,d + 'np) 

G 


( 9 ) 


which represents the relaxation time of a droplet with charac¬ 
teristic size H, determined by viscous and Capillary forces. Ei- 
nally, in figures I(e)-(f) we benchmark our numerical simula¬ 
tions under the effect of a steady elongational fiow, = £z, 
Ux = —ex/2, Uy = —Pyjl, with e the elongation rate. The as¬ 
sociated normal stress difference f{rp){Czz — Cxx) developed 
in steady conditions can also be evaluated analytically [39,42] 
and exact expressions are here omitted for simplicity. Similarly 
to the shear rheology, we have carried out numerical simula¬ 
tions in a three dimensional cubic domain of edge H consisting 
of H xH xH = 20x20x20 cells. Periodic conditions are now 
applied in all directions. The elongational rate is changed in the 
range 10“^ < £ < 10“^ Ibu and the polymer relaxation time in 
the range 10^ < Tp <10^ Ibu. Again, the various quantities are 
made dimensionless with the viscosity rfp and the relaxation 
time Tp, and they are plotted as a function of the dimension¬ 
less extensional rate A^ = Tp£ for different values of L^\ for 
small Ae the elongational viscosity is three times the polymer 
viscosity, while at large A^ we approach another constant value 
dependent on the finite extensibility parameter L^. 

Numerical simulations for the flow-focusing geometry will be 
performed with = 10"^: we note that the value of chosen 
rules out important thinning effects up to A ^ yTp ^ I — 10, 
which are above the values achieved in our simulations. We 
emphasize that the choice of a finite (but large) is the one 
that maximizes the effects of viscoelasticity, while keeping the 
computation stable. Such choice is indeed instrumental to en¬ 
sure upper bounds for the polymer elongational viscosity in 
steady elongational fiows, that would tend to diverge for the 
Oldroyd-B limit as soon as dimensionless elongational rates of 
order one are used (see figures I(e)-(f)). Obviously, inspecting 
the importance of by performing simulations of the flow- 
focusing geometry with other values of [15,18] is surely 
worth of investigation, but outside the scope of the present pa¬ 
per. 
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(c) Polymer first normal stress difference in shear flow 



(e) Polymer normal stress difference in elongational flow 




A 


(d) Polymer first normal stress coefficient in shear flow 



(f) Polymer elongational viscosity 


Fig. 1: Polymer shear and elongational rheology. Panels (a)-(b): we plot the shear stress and the shear viscosity (both scaled 
to the polymer viscosity r(p and polymer relaxation time Tp, see (l)-( 2 ) and text for details) as a function of the dimensionless 
shear A = T^yin a steady shear flow with intensity 7 . Symbols are the results of the numerical simulations [40,41] with different 
imposed shears, different Tp, and different values of L^. Numerical simulations for the flow-focusing geometry will be performed 
with = 10 "^ (see text for details), while the other values of are here used to give an idea of how much affects the 
viscoelastic response. The solid lines are drawn from theoretical predictions based on equations ( 6 ) and (7). Panels (c)-(d): we 
plot the dimensionless first normal stress difference and the first normal stress coefficient. Panels (e)-(f): we plot the dimensionless 
normal stress difference and elongational viscosity developed in steady elongational flow (see text for details) as a function of 
the dimensionless elongational rate A^ = Tp£. In all cases we also report the prediction for Hookean dumbbells (Oldroyd-B limit, 
> 1 ). 
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Ca 

Q 

Lx xLyXH 
cells 

Id 

Ibu 

he 

Ibu 

riP 

Ibu 

Tp 

Ibu 

De 

0.0056 

0.5-4 

1150x250x50 

1.75 

1.75 

0.00 



0.0056 

0.5-4 

1150x250x50 

1.05 

1.75 

0.69 

5 - 80 X 10^ 

0.2-3.2 

0.0056 

0.5-4 

1150x250x50 

1.75 

1.05 

0.0 

5 - 80 X 10^ 

0.2-3.2 

0.007 

0.5-4 

1150x250x50 

1.75 

1.75 

0.00 



0.007 

0.5-4 

1150x250x50 

1.05 

1.75 

0.69 

5-80x10^ 

0.2-3.2 

0.007 

0.5-4 

1150x250x50 

1.75 

1.05 

0.0 

5 - 80 X 10^ 

0.2-3.2 


Table 1: Parameters for the numerical simulations with flow-focusing geometry: Ca is the Capillary number, which quantifies the importance 
of the viscous forces with respect to the surface tension forces, and Q is the fiow-rate ratio of the two fluids. The fiow-focusing geometry 
is embedded in a rectangular parallelepiped with size Lx'xLyX H, where H is the characteristic edge of the channels with square {H x H) 
cross section (see also figure 2 for a sketch). The dynamic viscosity of the Newtonian solvent inside the dispersed phase is T]^, while ric is the 
dynamic viscosity of the Newtonian solvent inside the continuous phase, rjp is the viscosity of the polymers in the dispersed or continuous 
phase. The polymer relaxation time is Tp, while De is the Deborah number based on definition (8). 



4 Results and Discussions 

In this section we report the results for the break-up of liq¬ 
uid threads in the confined flow-focusing geometry at fixed 
Capillary number. We aim at exploring the effects of both 
MV and DV on the scenarios which are known for Newtonian 
phases. In particular, Liu & Zhang [10] performed LB simu¬ 
lations and systematically characterized three typical flow pat¬ 
terns for Newtonian phases, dependently on the flow-rate ratio 
Q. These flow patterns are actually well reproduced by our nu¬ 
merical simulations, as shown in flgure 3 for Ca = 0.007. At 
very low flow-rate ratio Q, the droplets are formed at the cross¬ 
junction (DCJ) due to the squeezing mechanism [7,47]. Panels 
(a)-(d) show the liquid thread prior to the first, second, third 
and fourth break-up at time t = to1 .6rehear, t = to2 .3T,hear, 
t = to -h 3.0T,hear, ^ufl ^ + 3.8 T, hear, respectively. Notice that 

we have used the characteristic shear time T,hear = H/vc as a 
unit of time, while is a reference time (the same for all 
simulations). As we can see, the action of the co-flowing liq¬ 
uid periodically breaks droplets and the break-up point is sta¬ 
ble in time. Upon increasing Q, droplets are found to pinch- 
off downstream of the cross-junction (DC). Panels (e)-(h) re¬ 
fer to a flow-rate ratio 2 = 3.0 and show the liquid thread 
at time t = to ^ T It, hear, t = to ^ 1.6T,hear, t = to ^ 2. It, hear, and 
t = to 2.5T,hear, respectively. It is evident that the dispersed 
thread actually becomes unstable after a distance from the 


cross-junction and this distance increases as a function of time, 
which is a distinctive signature of the DC regime. Practically, 
Lij is computed as the distance between the break-up point and 
the comer up-stream at the cross-junction. As the flow-rate ra¬ 
tio Q increases to a critical value, stable parallel flows (PF) are 
observed, where the three incoming streams co-flow in parallel 
downstream of the cross junction without pinching. This is ev¬ 
idenced by Panels (i)-(l), showing the liquid thread dynamics 
at time t = to-h 0.6T,hear, t = to ^ 1.2T,hear, ^ = ^0 + 1.8T,hear, and 
t = to 2.4T,hear, respectively. In this latter case, the flow-rate 
ratio is 2 = 4.0. In addition, the transitions from DCJ to DC 
and from DC to PF are influenced by the Capillary number. As 
Ca increases, the threshold value of flow-rate ratio at which the 
transition occurs decreases, and the width of the DC regime 
also decreases [10]. 

4.1 Effects of Matrix Viscoelasticity (MV) 

To appreciate the effects of matrix viscoelasticity in the DCJ 
regime, in flgure 4 we show the break-up process for different 
De. We report the density contours of the dispersed phase for 
three cases: (a) Newtonian matrix (no feedback stresses) at time 
t = to + 3.0T,hear; (A) sflghtly viscoelastic matrix with De = 0.2 at 
time t = ^0 + 2.8T,hear; (c) viscoelastic matrix with Deborah num¬ 
ber just above unity (De = 2.0) at time t = to 2.8T,hear- All the 
other flow parameters are kept flxed, Ca = 0.007, Re = 0.018, 
A = 1 and Q= 1.0. The mechanism of break-up is essentially 
the same, since droplets are formed at the cross-junction due to 
the squeezing mechanism. We observe, however, that as soon 
as the degree of viscoelasticity sensibly increases (De = 2.0), 
droplets pinch off with a slightly smaller size (flgure 4(c)) 
and the shape of the thread prior to break-up is different from 
that of the Newtonian case (flgure 4(a)). The opportunity of¬ 
fered by numerical studies is to allow for a systematic anal¬ 
ysis of the various forces present in the equations of motion 
(see Eqs. (1) and (2)). We can use this advantage to go deeper 
into the problem of the break-up in presence of viscoelastic¬ 
ity and visualize the magnitude of the polymer feedback stress 
i^f{rp)^ in Eq. (1)) for the two MV cases presented in flg¬ 
ure 4 at time t = to 2.8T,hear- In flgure 5 we report the density 
contours of the dispersed phase overlaid on snapshots of the 
magnitude of the polymer feedback stress. Working with MV, 
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Fig. 3: Dynamics and Break-up of liquid threads in the flow-focusing geometry from simulations with Newtonian phases. The 
dispersed (d) phase (with viscosity r]j) enters at the inlet of the main channel with average velocity vj and periodically breaks due 
to forces created by the cross-flowing continuous (c) phase from the side channels, the latter characterized by a viscosity and 
average velocity v^. We consider the simplest case of square ducts with edge H. The various hydrodynamical parameters can be 
grouped into the following dimensionless numbers: the Capillary number calculated for the continuous phase, Ca = <7, the 

Reynolds number Re = pVcHjric^ the viscosity ratio A = PdMe^ and the flow-rate ratio Q = Vdlvc = Qd!Qc^ where Qd = VdH^ 
and Qc = VcH^ are the flow-rates at the two inlets. The case reported in this flgure corresponds to Ca = 0.007, Re = 0.018, 
A = 1. The first row shows results for 2 = 1.0 (DCJ regime); the second row corresponds to 2 = 3.0 (DC regime); the third row 
corresponds to 2 = 4.0 (PF regime). The various regimes (DCJ, DC, PF) are described in the text. In all cases we have used the 
characteristic shear time T,hear = as a unit of time, while is a reference time (the same for all simulations). 





Fig. 4: Effects of matrix viscoelasticity (MV) on droplet formation in the flow-focusing geometry at fixed Ca = 0.007, Re = 
0.018, A = 1 and 2=1-0 (DCJ regime, see text for details). Three cases are compared: (a) Newtonian case (De = 0.0) at time 
r = ^0 + 3.0T,hear; (b) sUghtly MV case with De = 0.2 at time r = + 2.8T,hear; (c) MV case with De = 2.0 at time r = + 2.8T,hear- 

In all cases we have used the characteristic shear time T,hear = as a unit of time, while is a reference time (the same for all 
simulations). 


the feedback stress is obviously non zero only in the contin- the fluid, are well pronounced in correspondence of the cor- 
uous phase. These images clearly show that the stretching ef- ners where the side channels meet the main channel. This is a 
feet on the polymers, and consequently the feedback stress on consequence of the flow structure and the associated stream- 
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(a) t = to + 3.0T,hear, De = 0.0 (b) t = to + 2.8^, De = 0.2 (c) t = to + 2.8^, De = 2.0 


Fig. 5: Effects of matrix viscoelasticity (MV) in the DCJ regime. Density contours of the dispersed phase are overlaid on the 
polymer feedback stress magnitude (see Eq. (1)) for different values of the Deborah number De at fixed Ca = 0.007, Re = 0.018, 
A = 1 and 2 = 1.0. The value of the fiow-rate ratio Q is such that the squeezing dominated-DCJ regime [10] drives the break-up 
process. By definition, the Newtonian case (Panel (a)) has no feedback stress. As De is increased, we see that the fiow in the 
continuous phase develops enhanced polymer feedback stress in the comers where the side channels meet the main channel. 
The feedback is higher downstream at the cross-junction, causing enhanced droplet deformation and quantitatively changing the 
break-up process (see also figure 8). In all cases we have used the characteristic shear time T,hear = as a unit of time, while 
to is a reference time (the same for all simulations). 


Force Effective X Force Effective Y 



(a) t = to + 2 . 8 Tshear, De = 2.0 (b) t = to+ 2 . 8 ^, De = 2.0 


Fig. 6: Panels (a)-(b): x and y component of the effective force (see Eq. (10)) for a matrix viscoelasticity (MV) case at 
t = to + 2.8T,hear, De = 2.0, Ca = 0.007, A = 1 and Q= 1.0. When the liquid thread obstructs the main channel, we observe 
the formation of an elastic boundary layer around the emerging thread. This promotes an extra force which points towards the 
comers downstream of the emerging thread and triggers an extra deformation at the cross-junction. In all cases we have used the 
characteristic shear time T,hear = as a unit of time, while ^o is a reference time (the same for all simulations). A vectorial plot 
of the effective force is provided in figure 7. 


lines which provide stretching of the polymers in those re¬ 
gions. As De is increased, we see that the fiow in the contin¬ 
uous phase develops enhanced polymer feedback stress. The 
enhancement of the polymer feedback stress has quantitative 
effects on the velocity field, as predicted by equation (1). One 
has also to remark that viscoelastic forces provide a contribu¬ 
tion to the shear forces. This happens in simple shear fiows 
and also for weak viscoelasticity [42,43], where we expect that 
the viscoelastic stresses closely follow the viscous stresses, i.e. 
• [f{rp)^] K, V • {rip{Vuc + (Vmc)^))- Obviously, this 


cannot be the case when viscoelasticity is enhanced and the 
Deborah number is above unity. For this reason, to better vi¬ 
sualize the importance of the viscoelastic forces in comparison 
with the Newtonian case, one can define an effective force (F^ff) 
as [41] 

^ V ■ [/(rp)^] - V (r/p( Vuc + (Vn,)^)). (10) 

Tp 

Since all our simulations are performed with the same shear 
viscosity, the effective force gives an idea of how much the 
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Fig. 7: A vectorial plot of the effective force given in figure 6 
is provided. 



Fig. 8: Effects of matrix viscoelasticity (MV) on the droplet 
length. We report the droplet length (Lj) normalized to the 
characteristic edge of the channels (H) versus the fiow-rate ra¬ 
tio Q for different values of De for fixed Ca = 0.007, Re = 
0.018, and A = 1. We see that the droplet length is linearly in¬ 
creasing with fiow-rate ratio Q for all the three cases presented 
in the figure (Newtonian, slightly viscoelastic (De = 0.2) and 
viscoelastic (De = 2.0) case). Half of the interface thickness 
has been used as errorbar. Linear fits are drawn from the ex¬ 
pected behaviour (11), with ai^2 two fitting parameters . 


viscoelastic system differs from the corresponding Newtonian 
system with the same viscosity. If present (F,ff ^ 0), this change 
is solely attributed to viscoelasticity. For the case reported in 
figure 5(c), we have analyzed the effective force in the xy-plane 
at z = i7/2. Results are reported in figure 6, where we show two 
distinct plots for the x and y component of F,ff. The figure re¬ 
veals the formation of an elastic boundary layer and the emer¬ 
gence of extra forces which point towards the comers down¬ 
stream of the emerging thread. A vectorial plot of the effective 
force is provided in figure 7 to highlight how the direction of 
the force varies across the junction. 

To proceed further and complement the results of figures 4 
and 5, we have characterized the effect of MV and the fiow- 


rate ratio on the droplet length Lj. In figure 8 we study Lj 
normalized to the characteristic edge of the channels (H) for 
fixed Ca = 0.007, Re = 0.018, and A = 1. We report L^/H as 
a function of Q, for various values of the Deborah number. In 
the squeezing-dominated DCJ regime, a linear relation between 
the droplet length Lj and fiow-rate ratio Q is expected [47] 


— = {ai + a2Q) ( 11 ) 

where ai^2 are coefficients of order one. The specific values of 
these coefficients depend on the Capillary number and channel 
geometry [2,10,48,49]. The study that is probably best useful 
for us is that of Liu & Zhang [10], who computed the coeffi¬ 
cients a\^2 for a wide range of channel geometries and Cap¬ 
illary numbers in the Newtonian case. For a square duct with 
Ca = 0.007, their prediction (p. 11 in [10]) is ai = 1.17 and 
a 2 = 0.85. From a best fit of our data, we obtain a\ = 1.035 
and a 2 = 0.48. The value of a\ is actually in the ballpark of 
that found by Liu & Zhang [10], whereas our a 2 is smaller. 
As for the value of a 2 , we notice that Liu & Zhang define Q 
with half of the average velocity in one of the two side chan¬ 
nels (see figure 3 in [10]), which can explain why our fitted 
value of a 2 is roughly a factor 2 smaller than what they pro¬ 
vide. Other discrepancies may be attributed to different factors, 
including a different Reynolds number Re [50] and viscosity 
ratio A, although we do not expect a major infiuence of the lat¬ 
ter in the squeezing regime. Moreover, the resolution used by 
Liu & Zhang [10] is smaller than ours, although the authors 
acknowledge a grid independence with several different fiow 
conditions and find that the mesh refinement produces varia¬ 
tions not more than 5%. 

For the viscoelastic cases, the relation between the droplet 
length Ld and fiow-rate ratio Q is still linear, but an overall 
decrease is observed when comparing with the corresponding 
Newtonian data. In particular, the coefficient a 2 is not sensi¬ 
bly different, whereas we acknowledge a decrease in a\. We 
go from a\ = 1.035 in the Newtonian case to a\ = 0.8 when 
De = 2.0. This behaviour can be better understood if we ex¬ 
plore the basic ingredients which lie at the core of the observed 
linear scaling law [47]. The break-up is the result of two dis¬ 
tinct physical processes: first, a blocking process, where the 
thread of the dispersed phase grows until it effectively blocks 
the cross-section of the side channel. At this particular mo¬ 
ment, the blocking length Lbiock of the plug is of the order of 
the channel width, say aiH (with ai a constant of order unity). 
Afterwards, the necking process starts. For a neck with a char¬ 
acteristic width a 2 H (a 2 is a constant, again, of order unity) 
and squeezing at a rate approximately equal to the average ve¬ 
locity (Qc/H^), it takes a time T ^ a 2 HH^/Qc to complete 
the squeezing process. During this time, the plug continues to 
elongate with a rate Qd/H^. The resulting “squeezing length” 
of the plug is therefore ~ = ^liiQdlQc^ Con¬ 

sequently, the final dimensionless length Ld/H of the droplet 
can be expressed as ^ ai + OC 2 Q. It is clear from fig¬ 

ure 8 that viscoelasticity directly affects the offset constant a \, 
i.e. it changes the blocking process of the squeezing regime. 
This is due to the extra force generated by the polymers feed¬ 
back stresses which is visible in figure 6. This force effectively 
pulls the droplet interface downstream at the cross junction, so 
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(c) t = to + 2.9 w, De = 2.0 



Fig. 9: Effects of matrix viscoelasticity (MV) in the DC regime. Droplet formation in the flow-focusing geometry at flxed 
Ca = 0.007, Re = 0.018, A = 1 and Q = 3.0. Three cases are compared: (a) Newtonian case (De = 0.0) at time t = to + 3.0T,hear; 
(b) slightly viscoelastic case with De = 0.2 at time t = to + 2.8T,hear; (c) MV case with De = 2.0 at time t = to + 2. 9 T,hear- In all 
cases we have used the characteristic shear time T,hear = Hjvc as a unit of time, while to is a reference time (the same for ah 
simulations). 



{ 


(a) t = to + 2.8T,hear, Dc = 0.2 


(b) t = to + 2.9T,hear, Dc = 2.0 


Fig. 10: Effects of matrix viscoelasticity (MV) in the DC regime. Density contours of the dispersed phase overlaid on the polymer 
feedback stress magnitude (see Eq. (1)) for two characteristic values of the Deborah number De at flxed Ca = 0.007, Re = 0.018, 
A = 1 and 2 = 3.0. As De is increased, we see that the flow in the continuous phase develops enhanced polymer feedback stress 
at the junction between the side channel and the main channel. The feedback is higher downstream at the cross-junction. The 
effect of viscoelasticity is to promote a shift of the break-up point from downstream of the cross-junction to the cross junction 
itself. In ah cases we have used the characteristic shear time T,hear = F/^/vc as a unit of time, while to is a reference time (the same 
for ah simulations). 


that the blocking length, i.e. the distance that the droplet has to 
travel to obstruct the side channels, is effectively smaller. 

We next go on by quantifying the effects of polymers in the 
DC regime. As discussed in flgure 3, a distinctive feature of 
the DC regime, is that the break-up point moves progressively 
into the main channel where the dispersed phase breaks-up af¬ 
ter traveling a distance from the cross-junction. The break¬ 
up point is not stable, i.e. the distance increases as a function 
of time (see flgure 3(e)-(h)). In flgure 9 we can appreciate the 
effects of MV in the DC regime. In particular, we report den¬ 
sity contours of the dispersed phase for three cases: (a) New¬ 
tonian matrix at time t = to (b) slightly viscoelastic 
matrix with De = 0.2 at time t = to ^ 2.8T,hear; (c) viscoelas¬ 
tic matrix with Deborah number just above unity (De = 2.0) 
at time t = to 2.9T,hear- Ah the other flow parameters are kept 
flxed, Ca = 0.007, Re = 0.018, A = 1 and Q = 3.0. The time 
of the observation is essentially the same (« + 2.8T,hear) and 
the effect of viscoelasticity is to promote a shift of the break¬ 
up point from downstream of the cross-junction (flgure 9(a)) 


to the cross junction itself (flgure 9(c)). The density contours 
of the dispersed phase overlaid on snapshots of the magnitude 
of the polymer feedback stress are reported in flgure 10. Again, 
we attribute the observed behaviour to the extra force generated 
by the polymer feedback stress at the cross-junction: this is ef¬ 
fectively perturbing the droplet shape, provoking the break-up 
immediately after the cross junction. In some sense, this effect 
may be translated in an enhanced stability of the DCJ regime, 
a fact that will be better quantifled in flgure 13. 

Finally, we investigate the effect of MV on the PF regime. This 
is shown in flgures 11 and 12, where we can appreciate how the 
perturbation induced by viscoelasticity to the droplet shapes at 
the cross-junction are enough to spoil the PF itself leading to 
droplet break-up (see flgure 11(c)). 

To provide an overview of the results for the various regimes, in 
flgure 13 we report the first break-up distance normalized to the 
characteristic edge of the channel (Lh/H) as a function of the 
flow-rate ratio Q for various Deborah numbers. The other pa¬ 
rameters are kept flxed to Ca = 0.007, Re = 0.018, and A = 1. 
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Fig. 11: Effects of matrix viscoelasticity (MV) in the PF regime at fixed Ca = 0.007, Re = 0.018, A = 1 and Q = 4.0. Three 
cases are compared: (a) Newtonian case (De = 0.0) at time r + 2.4T,hear; (b) slightly viscoelastic case with De = 0.8 at time 
t = to ^ 2.4T,hear; (c) viscoelastic case with De = 2.0 at time ^ + 2.4T,hear- In all cases we have used the characteristic shear time 

'Z^shear = H jvc di Unit of time, while is a reference time (the same for all simulations). 



(a) t = to + 2.4T,hear, De = 0.8 (b) t = to+ 2.4^, De = 2.0 


Fig. 12: Effects of matrix viscoelasticity (MV) in the PF regime. Density contours of the dispersed phase overlaid on the polymer 
feedback stress magnitude (see Eq. (1)) for two characteristic values of the Deborah number De at fixed Ca = 0.007, Re = 0.018, 
A = 1 and Q = 4.0. Similarly to figure 10, as De is increased, the fiow in the continuous phase develops enhanced polymer 
feedback stress at the junction between the side channel and the main channel. This results in the appearance of a break-up point 
downstream at the cross-junction. In all cases we have used the characteristic shear time T,hear = as a unit of time, while 
is a reference time (the same for all simulations). 


In the Newtonian case, the DCJ regime {Lb/H ^ 1) is valid for 
small 2 up to 2 ^ 2.0 [9,10], whereas at larger Q we observe 
the emergence of the DC regime, with the break-up distance 
that is not stably located at the cross-junction anymore (circles 
indicate unstable L^), but moves progressively downstream of 
the cross-junction. For the viscoelastic cases, the emergence 
of the DC regime is delayed. Just to give some numbers, an 
increase of De from De = 0.0 to De = 2.0, results in the emer¬ 
gence of the DC regime for Q ^ 3.0 instead of 2 ~ 2.0. To 
stress even more the stabilizing character of the MV in the 
transition from DCJ to DC, at fixed fiow-rate ratio 2 = 2.0, 
we report in figure 14(a) the break-up distances up to the fifth 
break-up as a function of De. A clear distinction emerges by 
comparing the situation at De = 0 (Newtonian case) and that 
with De just above unity, where the break-up point is mani¬ 
festly stabilized. For fixed Q = 3.0, we report in figure 14(b) 
the first break-up distance as a function of De, showing a linear 
decrease, = L^(De = 0) — aDe, with a a constant of order 
unity. 


4.2 Effects of Droplet Viscoelasticity (DV) 

In this subsection we will discuss the effect of DV for the 
break-up of threads in the confined flow-focusing geometry. 
We will be mainly addressing the role of DV in the DCJ regime, 
where quantitative results for the MV have been obtained as a 
function of the fiow-rate ratio Q (see figure 8). In figure 15(a) 
we study droplet length Lj normalized to the characteristic 
edge of the channels (H) for fixed Ca = 0.007, Re = 0.018, 
and A = 1. A linear relation between droplet length Lj and 
fiow-rate ratio Q is again observed, which echoes the results 
for MV cases in figure 8. For a given fiow-rate ratio Q, initially 
the droplet length Lj decreases with De but, when De is above 
unity (De = 2), we find that Lj is slightly larger compared to the 
Newtonian case. These changes are however small if compared 
to those achieved for the MV cases. To quantitatively compare 
the effect of DV and MV we have collected in figure 15(b) data 
from figures 8 and 15(a). In figure 16 we compare both DV and 
MV, by studying normalized to the characteristic edge of the 
channels (H) for fixed Ca = 0.007, Re = 0.018, and A = 1. In 
both cases, the effect of viscoelasticity is a stabilizing one but, 
again, it is more pronounced in the case of MV. 
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Fig. 13: Effects of matrix viscoelasticity (MV) on the break¬ 
up point. We report the break-up distance of the first break¬ 
up (Lij) normalized to the characteristic edge of the channels 
(//) for fixed Ca = 0.007, Re = 0.018, and A = 1 for different 
De and Q. As we can see, as soon as De is just above unity, 
the fiow undergoes a transition and the break-up distance 
moves towards the cross-junction. Black open circles represent 
the cases where the break-up length is changing with time and 
moves progressively downstream of the cross-junction. 


for a fixed L^, i.e. for a fixed maximum elongation of the poly¬ 
mers. The parameter is actually related to the extensional 
viscosity of the polymers [39,42]. Similarly to the work that 
we recently developed for characterizing droplet deformation 
and break-up in confined shear fiow [41], it would be impor¬ 
tant to study the impact of a change in the finite extensibility 
parameter for the geometries studied in this paper [15,18]. 
Finally, we wish to underscore the role played by numerical 
simulations: the numerical simulations are crucial to elucidate 
the relative importance of the free parameters in the model, and 
to visualize the distribution of the polymer feedback stresses 
during the motion of the droplets. Thanks to these insights, it 
was possible to correlate the distribution of the stresses to the 
droplet shape and the corresponding break-up morphology. 
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Framework Programme (FP7/2007-2013) / ERC Grant Agree¬ 
ment N. 279004. We acknowledge support from the COST- 
Action MP1305 and the computing hours from ISCRA B 
project (POLYDROP), CINECA Italy. MS acknowledges use¬ 
ful discussions and fruitful exchange of ideas with Prof. H. Liu 
during his visit in June 2014. 


5 Conclusions 


6 Hybrid Lattice Boltzmann (LB) Models - 
Finite Difference (FD) Schemes 


A crucial point for designing, developing and exploiting micro- 
and nanofiuidic devices is to achieve a systematic control over 
the process of formation of droplets as a function of the 
materials and fiow parameters. The confinement that naturally 
accompanies fiows in small devices has significant qualitative 
and quantitative effects on the droplet dynamics and break-up. 
Moreover, relevant constituents have commonly a viscoelastic 
- rather than Newtonian - nature. One of the most common 
droplet generator is the flow-focusing geometry, based on 
the focusing of a stream of one liquid in another co-flowing 
immiscible liquid. In this paper we have presented results 
based on three dimensional mesoscale numerical simulations 
to highlight the non trivial role played by confinement and 
non-Newtonian effects. We have worked in the flow conditions 
previously studied by other authors [9,10,51], where New¬ 
tonian droplets in a Newtonian matrix are squeezed either at 
the cross-junction, or form downstream of the cross-junction 
(DC). The transition between these two regimes has been 
found to be affected by viscoelasticity, and in particular by 
those situations where viscoelastic properties are confined in 
the continuous phase. This is due to the fact that the action 
of the co-flowing liquid that periodically breaks droplets is 
supplemented with the polymer feedback stresses that are 
well pronounced in correspondence of the comers where the 
side channels meet the main channel. The enhancement of 
the polymer feedback stress promotes an extra force which 
points towards the comers and inevitably perturbs the flow 
properties. Specifically, an extra deformation of the droplet 
at the cross-junction is introduced and is able to trigger the 
break-up immediately after the cross junction itself, thus 
resulting in a delayed DCJ to DC transition. 

We remark that the results presented in this paper are obtained 


In this appendix we report some details of the numerical 
scheme used. The interested reader can find more extensive 
technical details in a dedicated paper [41]. We use a hybrid al¬ 
gorithm combining a multicomponent Lattice-Boltzmann (LB) 
model with Finite Differences (FD) schemes. LB is used to 
model the macroscopic hydrodynamic equations, while FD is 
used to model the dynamics of the polymers. The LB time dy¬ 
namics considers the discretized probability density function 
to find at position r and time t a fluid particle of com¬ 
ponent a =A^B with velocity c/. The dispersed (d) and the con¬ 
tinuous (c) Newtonian phases in Eqs. (1) and (3) are character¬ 
ized by a majority of one of the two components, i.e. majority 
of A (B) in the dispersed (continuous) phase. The LB evolution 
scheme with a unitary time-step reads as follows [21]: 

fai{r + Ci,t + l)-fai{r,t) = Y^^ij{faj-/af) + A^„j. ( 12 ) 

j 

The first term in the rhs of Eq. (12) is the (linear) collisional op¬ 
erator, expressing the relaxation of fat towards the local equi¬ 
librium f^f^. All the numerical simulations make use of the 
D3Q19 model (3d with 19 velocities) whose discrete velocity 
set is given by 


{ ( 0 , 0 , 0 ) / = 0 

(± 1 , 0 , 0 ),( 0 ,± 1 , 0 ),( 0 , 0 ,± 1 ) /= 1...6 . 

(±1,±1,0),(±1,0,±1),(0,±1,±1) / = 7...18 

(13) 

The expression for the equilibrium distribution is [52,53] 


Aeq) _ ^ 

Jai — ^iPa 


UCi 


uu : (ciCi-l) 


2cf 


(14) 
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Fig. 14: Effects of matrix viscoelasticity (MV) on the break-up distance. Panel (a): we report the break-up distances up to the 
fifth break-up as a function of De for fixed Q = 2.0, Ca = 0.007, Re = 0.018, and A = 1. The break-up distance is normalized 
to the characteristic edge of the channels (H). Panel (b): we show the break-up distance of the first break-up as a function of the 
Deborah Number for a fixed flow-rate ratio 2 = 3.0. 




(a) Ca = 0.007 


(b) Ca = 0.007 


Fig. 15: Effects of matrix viscoelasticity (MV) and droplet viscoelasticity (DV) on the droplet length. Panel (a): we report the 
droplet length (Lj) normalized to the characteristic edge of the channels (H) versus the fiow-rate ratio Q for different values of 
De in the DV case. All the other parameters are kept fixed to Ca = 0.007, Re = 0.018, and A = 1. Panel (b): both the cases with 
MV and DV are compared. Half of the interface thickness has been used as errorbar. Linear fits are drawn from the expected 
behaviour (11), with ai^2 used as fitting parameters. 


and the weights Wi are 



1/3 

1/18 

1/36 


/ = 0 

/ = 1 ... 6 , 
/ = 7...18 


(15) 


where Cs is the isothermal speed of sound and u is the fiuid ve¬ 
locity. The operator in Eq. (12) is the same for both com¬ 
ponents and has a diagonal representation in the mode space: 
the basis vectors H]^ (k = 0,..., 18) of the mode space are used 
to calculate a complete set of moments, the so-called “modes” 
Mak = HiHkifai (k = 0,..., 18) [52-55]. Hydrodynamic vari¬ 


ables are obtained from the lowest order modes: the density of 
both components and the total density are pa = ^aO = Hifai^ 
P = = JlaPa, while the next three moments = 

(mai^mai^fnas), define the velocity of the mixture 



a 


2p 


1 

P 


a i 


2p‘ 


(16) 


The higher order modes refer to the viscous stress tensor, 
and also other modes (the so called “ghost modes”) which 
do not appear in the hydrodynamic equations. Since the op¬ 
erator ^ij possesses a diagonal representation in mode space, 
the collisional term describes a linear relaxation of the modes, 
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In the above equations, p = Y.aPa = JLa^^Pct represents the 
internal (ideal) pressure of the mixture. The quantity Da is the 
diffusion flux 


Da — P 



with p the mobility coefficient: 




( 22 ) 


(23) 


As for the internal forces, we will use the “Shan-Chen” 
model [44,56,57] for multicomponent fluids 


Fig. 16: We report the break-up distance of the first break¬ 
up (Lh) normalized to the characteristic edge of the channels 
(H) versus the flow-rate ratio Q for different values of De. Re¬ 
sults for droplet viscoelasticity (DV) are here compared with 
those of matrix viscoelasticity (MV) already reported in figure 
13. Black open circles represent the cases where the break-up 
length is changing with time and moves progressively down¬ 
stream of the cross-junction. All the cases reported have fixed 
Ca = 0.007, Re = 0.018, and A = 1. 


= {l-\-Xk)mak-^J^ak^ where the “post” indicates the post- 
collisional mode and where the relaxation frequencies are 
related to the transport coefficients of the modes (see also be¬ 
low). The term is the k-th moment of the forcing source 
A^f. such term embeds the effects of a forcing term with den¬ 
sity Qa [52,54]. The term g = Y.a 9a in Eq. (16) relates to the 
total force. The forcing source is tuned in such a way that the 
correct hydrodynamic equations are obtained [55]: given the 
relaxation frequencies of the momentum (—Am), bulk (—A/,) 
and shear (—A^) modes, the forcing source reads 


Tg: (c;C;-C?l) , 

(17) 

(18) 

LB reproduces the continuity equations and the NS equations 
for the total momentum [52,54,55] with shear viscosity Ps and 
bulk viscosity r/^: 


^ 


Wi / 2 -j- Am \ 


V 2 


Wi 

9a • Q + ^ 


+ V • (Pau) — V • Dct: (19) 


p [dtU-\- {u-'V)u] = — V/7 


Ps ( V'U + (V'U)^ - -I(V-'U) ) -^Pbt{V-u) 


( 20 ) 


The relaxation frequencies of the bulk and shear modes in (12) 
are related to the viscosity coefficients as 


ns = -pc] 





( 21 ) 


Pair) = -gABpair)'^ ^ Wipai{r + Ci)ci a,a'=A,B 

i a'^a 

(24) 

where gAB is a coupling parameter that regulates the inter¬ 
actions between the two components. When gAB is chosen 
above a critical value, phase segregation occurs with the for¬ 
mation of stable interfaces with a positive surface tension. 
Due to the effect of interaction forces, the internal pressure is 
modified by the “interaction” pressure tensor [58], i.e. 
P = p\^ , with 

P^"'\r) = ^gABpA{r)Y,"^iPB{r + Ci)ciCi 

^ ‘ (25) 

+ ^gABPB{r)'^WiPA{r + Ci)CiCi. 


Wetting properties are introduced with a tuning of the density 
in contact with the wall [59,60]. The numerical simulations 
presented are carried out with gAB = 1.5 Ibu in (24), corre¬ 
sponding to a surface tension cr = 0.1 Ibu and associated bulk 
densities pA = 2.0 Ibu and pB =0.1 Ibu in the A-rich phase. 
The relaxation frequencies in (21) are set to Am = —1.0 Ibu 
and A^ = A^, thus reproducing the viscous stress tensor given in 
Eqs. (1) and (3)). The viscosity ratio of the LB fluid is changed 
by allowing A^ to depend on space 

+ 0 = ?]</(/+('/'))+nc{/-(</»)) (26) 


where (j) = (l){r) 
sen as 


^ (PaW+PbW) - functions /±(0) are cho- 



l±tanh(0/O ^ 


(27) 


which allows to recover the Newtonian part of the NS equa¬ 
tions reported in Eqs. (1) and (3) with shear viscosities Pd and 
Pc. The smoothing parameter ^ is chosen sufficiently small so 
as to match analytical predictions on droplet deformation in 
presence of viscoelastic stresses (see [41] for all details). 

As for the polymer evolution given in Eq. (2), we follow the 
two References [61,62] to solve the FENE-P equation. The 
polymer feedback stress is computed from the FENE-P evolu¬ 
tion equation and used to change the shear modes of the LB [41, 
52,53]. The feedback of the polymers is modulated [32] in 
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space with the functions /±(0) 
p [dtu^{u • V)n] = — VP 

+ [(^j/+(0) + 'nc/-(0))(Vn + (Vn)^)] ^28) 

+ ^V[/(rp)^/±(<^))]. 

Tp 

By using /- (0), we recover a case where the viscoelastic prop¬ 
erties are confined in the continuous (c) phase, while the use 
of the function /+(0) allows to recover a case where the vis¬ 
coelastic properties are confined in the dispersed (d) phase. 
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